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An exact, nonlocal, finite step-size algorithm for Monte Carlo simulation of theories with dynamical fermions 
is proposed. The algorithm is based on obtaining the new configuration U from the old one U by solving the 
equation M(U)r] = U)M(U)rj, where M is fermionic operator, r\ is random Gaussian vector, and ui is random 
real number close to unity. This algorithm can be used for acceleration of current simulations in theories with 
Grassmann variables. A first test was done for SU(3) QCD with purely fermionic term in the action. 



1. Introduction. 

From the earliest days of Monte Carlo simula- 
tions in lattice field theory fermionic fields have 
caused annoying difficulty, which stems from their 
being anticommuting variables. It is not immedi- 
ately straightforward to put fermions on a com- 
puter, which is expected to manipulate numbers. 
This problem is only an algorithmic one, since for 
the most actions in use one can eliminate fermions 
by an analytic integration. However, the result- 
ing expressions involve the determinants of very 
large matrices, making the numerical simulations 
extremely expensive. 

Over the years many interesting tricks have 
been developed to circumvent this problem. The 
most often used algorithm, Hybrid Monte Carlo 
(HMC) [Q, is now considered to be the standart 
simulation tool. In recent years it was success- 
fully challenged by the Multiboson method, which 
has many attractive features, like the possibility 
of simulation of an odd number of quark flavors. 
The most popular version of Multiboson method, 
proposed by M. Liischer H, has been intensively 
studied by many authors (see || and references 
therein/thereon), and after various improvements 
was claimed to be competitive or even slightly 
better than HMC (see for example The 
Multiboson algorithm, proposed by A. Slavnov 
H), is less known, and was tested in ||. Besides 
HMC and Multiboson method, I like to mention 
some interesting ideas, which may finally result in 



efficient algorithms for simulation of theories with 
fermions. They are based on the polymer J7J and 
Jordan- Wigner occupation number Q represen- 
tations of fcrmion determinant, the direct evalu- 
ation of Grassmann integrals the separation 
of low and high eigenmodes of the Dirac operator 
fic| l , and the direct simulation of loop expansion 
by means of using the stochastic estimators 

Despite the variety of different approaches to 
the problem, there is a common impression that 
all existing algorithms remain inefficient [p^[ . 
This is clearly seen if we compare the compu- 
tational cost for the theories with dynamical 
fermions from the one side, and purely bosonic 
theories from the other side. In QCD even for 
heavy quarks the simulations are orders of magni- 
tude slower than the ones in the quenched regime 
(the approximation where the quark determinant 
is set to be equal to unity) . Such a situation sug- 
gests that one should not stop the attempts to 
obtain relatively cheap fermion simulation algo- 
rithm. 

In this paper I propose a new computational 
strategy for treating dynamical fermions in Monte 
Carlo simulations. It has the virtues of exact- 
ness, nonlocality and finite step-size of update []. 

1 There is a great probability, that a successful fermionic 
algorithm must be nonlocal and have a finite step-size of 
update. Indeed, HMC implies global changes in configu- 
ration space, but with infinitesimal step-size. This drasti- 
cally slows down the simulations, becouse the high energy 
barries become nearly impassable. Multiboson algorithm 
has the advantage of finite step-size, but it updates the 
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Unlike the other exact algorithms, the CPU-time 
required for the update in this algorithm grows 
only lineary with the volume V of the system. 
This makes the new method to be rather attrac- 
tive for practical simulations, although there are 
still some problems to be solved. 

2. The algorithm 

Suppose that we aim at sampling the partition 
function of theory with two flavours of degenerate 
fermion fields: 

Z ferm = J Det (m*[U]M[U]) dU (1) 

where M is the discretized fermion operator act- 
ing at some vector space fi, U denotes the bosonic 
degrees of freedom coupled to fermion fields. For 
simplicity of formulae I postpone the inclusion of 
purely bosonic action Sb[U]. The sampling of full 
partition function 

S f „^ J e-^> Det (M< mm ) 4U (2) 

will be considered later. 

At the beginning let me introduce a simple aux- 
iliary algorithm, which will guide the later con- 
struction and help to understand the basis of new 
method. Starting from old configuration U, one 
executes the following instructions: 

Prescription A 

• generate random vector r\ £ O with Gaus- 
sian distribution: Pg{v] = Zg 1 e _ '' ) 

• propose new configuration U with symmet- 
ric probability 

P Q [U^U] =P q [U^U] (3) 

• accept U with the probability 

PacdmU^ U] =mm(l,el"l 2 -l^ lM "l 2 )(4) 

fields locally. The large autocorrelation times suggest that 
such local changes are not very 'physical' in this approach. 



Here Zq normalizes the distribution: 

Z G = J e- M2 dr]drf (5) 

In eq.(|]) and below I use the following short no- 
tations: M ee M[U]\ M ee M[U] . 
The full transition probability 

P[U->U] ee P [U^U] x 

Jdqdrf PqM PaccbT, U^U\ (6) 

satisfies the detailed balance condition with re- 
spect to the partition function ([!]): 

Det(M^M)P[U -^U} = Det(M^M)P[U -> U] (7) 

This can be easily seen by making the change of 
variables 

r/^Af" 1 ?? ; 77+ -> ^(M" 1 )* (8) 

in the expression (||). Therefore, if the proposal 
matrix @ ensures ergodicity, then Prescription 
A is a valid algorithm for sampling the partition 
function (|l|). 

However, this exact algorithm is expensive due 
to the necessity to invert the large matrix M in 
the expression (Q). The cost of inversion is pro- 
portional to the volume of the system V, so if one 
updates the fields U locally, the computational 
cost of updating the entire configuration grows 
as 0{V 2 ). One can be more clever and perform 
the global updates by fixing 4> = Mr], introducing 
fictitious momenta p conjugate to the fields U (or 
some functions of U) , and pursuing a discrete in- 
tegration of Hamilton evolution (this is the core 
of HMC algorithm, which in its traditional form 
has a volume dependence like U 5 / 4 ). Neverthe- 
less, the procedure remains expensive, becouse 
one should tend the size of molecular dynamics 
steps to zero for minimizing the integration er- 
rors and again invert the matrix M during the 
integration. 

Can one somehow modify the Prescription A to 
escape the inversion of matrix M, implementing 
the global update U — > U in configuration space? 
The answer is 'y es ' ■ 

The key idea is to introduce into proposal ma- 
trix d3l) the dependency on 77 to make the new 
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configuration U satisfy the approximate equality 
M Mr) ~ rj, or equivalently, 



Mr/ ps Mr? 



(9) 



The symmetry of expression ([)]) under the inter- 
change of variables U <-> J7 will help to garantee 
the reversibility of algorithm. Moreover, the ac- 
ceptance probability (0) should be large, if the 
approximate equation (g) is close enough to its 
exact analogue. 

After these handwaving speculations let me 
present a rigorous construction. I do it in two 
steps. Firstly, I prove the detailed balance con- 
dition for some general Prescription B, in which 
the proposal matrix Po[u>, T] ; U— > U] depends on 
rj and random real number ui distributed in the 
narrow interval near unity. Secondly, I specify the 
particular choice of Pa[u),r) ; U—*U] and present 
the final algorithmic scheme. 

Prescription B 

• generate random vector r] £ f2 with Gaus- 
sian distribution: PgW\ — Zg 1 e~'' 7 ' 



• generate uj £ 



1 - e 



ability (J-[lu], satisfying 



1 

0J Z 



MM 



with the prob- 



(10) 



• propose new configuration U with the prob- 
ability Pq[u>, r\ ; U—>U], satisfying the fol- 
lowing symmetry relations: 



Po[u,ri; U^U] =P [-,v, U^U] (11) 



Po^M- 1 ^ ; U] = 
= P {- 1 M- 1 r r , U^U] 

UJ 



accept U with the probability 



(12) 



Here e is the algorithmic parameter lying in the 
interval < e < 1. 



The condition ( |10|) provides the invariance of 
the measure /i[uj]duj under the change of variable 
uj — > Indeed, using eq. (|To|) , one can easily 
check that for any integrable function f(u>) the 
following equality holds: 



fi[oj] /(w) dui 



li[u]f(l/u)du (13) 



Using condition (^lj) together with the prop- 
erty (|13|), one can demonstrate the symmetry of 
averaged over 77 , u> proposal matrix 



PS v [u 



U] 



fj,[uj]duj x 



dr)dtf P G [r)} P [uj,r) ; U U] 
under the interchange of variables U 

P£ v [u-*u] = ps v [u ->u] 



U: 



(14) 



(15) 



This makes the algorithm reversible. 

Finally, condition (|lj) ensures the fulfillment 
of detailed balance equation (^) for the full tran- 
sition probability 



P\U U] 



l-e 



n[uj]du> / dr/drf Pg[v] 



P a [uj, V ; U^U]P acc [ m U^U\ 



(16) 



Indeed, making the change of variables (||) in ex- 
pression (Il6|), one obtains: 



Bet (M f M) P[U -> U] 



PokM-'j,; U^U] (17) 

Expression ( p"7| ) is symmetric with respect to 
the interchange U <-> U due to the eqs.( |l2] . |l3"| ). 
Therefore, Prescription B is again a valid algo- 
rithm for sampling the partition function ([!]), if 
the averaged proposal matrix ( |l4| ) provides ergod- 
icity. 

Now I propose to specify the choice of the ma- 
trix Po[oj,r] ; {J— > U] by defining it through the 
equation: 



M 77 = u)Mr\ 



(18) 
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It means that, analytically or numerically, one 
finds some solution U of the equation ([l8]) and 
propose it as a new configuration. A good recipe 
for the numerical search may be the local iterative 
minimization of the quantity 

R = | (M - wM)r) | 2 (19) 

for fixed U,u>,r]. Starting from U = U, the mini- 
mization proceeds until R < S is reached, where 
S determines the accuracy of solving the eq.jl^). 
One should construct the minimization procedure 
in a way that ensures reversibility of the algo- 
rithm (i.e. probability to obtain the configuration 
U starting from U at any uo should be equal to 
probability to obtain U starting from U at 1/uj). 

It can be checked that the proposal matrix, de- 
fined ^ through the eq.(|l8|), satisfies the symme- 
try relation (fil|). Indeed, the equation (18) is 
invariant under the simultaneous interchange of 
variables U <-*•?/; u> «-► Therefore, the pro- 
posals Pq[ijJ, r)',U—> U] and Pq[1/u}, 77; C/ — > C/] are 
equiprobable, becouse they are defined through 
the same equation. 

The same logic is applicable for proving the 
fulfillment of symmetry relation (H). The lhs. 
of expression (|l^) is defined through the equa- 
tion MM~ l n = uji], meanwhile the rhs. is 
defined through the equation MM~ x n = jjT). 
These equations are both equivalent to the equa- 
tion M~ lr q = wM _1 ij. Therefore, the proposals 
P^uJ.M^n; U-> U] andP [l/w, M^m U -> U] 
are also equiprobable. 



Let us note, that on the surface (18) the accep- 
tance probability §Q) acquires the following simple 
form: 



Pa 



[ui, 77] = min (l, e^ 1 w 3 / ^ ) 



(20) 

The calculation of P acc becomes extremely cheap, 
since one does not need to invert the matrix M 
anymore. Moreover, the expression (|2(]) does not 
depend on U, so one can accept or reject u> (or 
the pair (u>,rj) ) even before solving the equation 

©• 

We also need to specify the probability fi[u>]. A 
good choice is the following expression: 

n[ui] oc mm(l, 1/w 2 ) (21) 

2 The definition means the choice of some concrete re- 
versible procedure of finding the solution U. 



which, evidently, satisfies the condition (|l0|). 
Considering together the expressions one 
gets the unified probability distribution for u: 



V[uj,rj\ = fi[u] x P acc [aj,n] 
oc min (l/w 2 , e^ 1- ^ 1 ) ) 



(22) 



Now we are ready to write down the final algo- 
rithmic scheme for sampling the partition func- 
tion (@): 

The algorithm 

• generate random vector rj £ fi with Gaus- 
sian distribution: Pg[v\ = Z^e"^ 



generate ui G 



1 - e 



1 

l-e 



with the prob- 



ability: V[u!,r]} oc min (l/w 2 , e*- 1 ^/' r ' ) 

• find the new configuration U by solving the 
equation: Mr] = u>Mtj 

The only computationally expensive ingredient 
of the algorithm is the obtaining of solution U of 
eq.([l8|). One can expect, that the usage of the lo- 
cal iterative minimization of the functional ( |l9| ) 
gives the computational cost proportional only 
lineary to the volume V. 

If the procedure for obtaining the solution of 
eq.([l8|) is fully specified, the algorithm has only 
one free parameter e, which controls the size of 
the deviation of u> from unity. A rather intriguing 
possibility is to set e = 0, therefore fixing 10 = 1. 
This makes the algorithm 'energy conserving', in 



a sense that |M _1 M?7| 2 



I '/I 



Unfortunately, 



in this case the search for the solution of eq.(J18|) 
by minimizing the functional (|l^) does not work, 
becouse one immediately obtains the trivial solu- 
tion U — U. However, in typical case (e.g. SU(3) 
QCD), the solutions of eq.(|l8|) are highly degen- 
erate, and the subspace of nontrivial solutions in 
configuration space is not empty. Finding the re- 
versible procedure of nontrivial solving of equa- 
tion Mn = Mn is the perspective subject for fu- 
ture investigations. 

In gauge theories one can use the gauge free- 
dom to simplify the procedure of solving the 
eq.([l8|). Indeed, under the gauge transformations 

(23) 
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eq.(|18|) acquires the form: 
Mrf = wMrf 



(24) 



where rj 9 = Grj, G is the matrix representing 
the gauge transformation on the vector space f2. 
Solving the eq. ( p4| ) can be particularly simple for 
some rf . 

Finally let me note, that the same algorithm 
can be used for simulations in theories with 
bosonic determinants, if we change the accep- 
tance (p(3). Sampling the partition function 



Det 



- 1 [m^[U]M[U] 



dU 



(25) 



one should use P acc [uj,r]] = rain (l,ev w J 
instead of the expression (20). 



3. Potential problems 

Despite the simple formulation and cheapness 
of the considered algorithm, it may be not ap- 
plicable for some models. The main danger is 
the possible absence of solutions of eq.(p^[). This 
problem may be principal (no solutions exist at 
all), or mild (no solutions exist for some partic- 
ular 77, u>). In the second case one can reject the 
pairs (77, uj) until the solution is found. The value 
of this 'hidden' acceptance would determine the 
actual efficiency of the algorithm in such a case. 

Another possible problem may be connected 
with the procedure of finding the solutions of 
eq. (|l8|) . Suppose that one uses the local iterative 
minimization of the quantity ( |l9|) for this pur- 
pose. If the functional ([l9|) can have some local 
minima at which R > 0, then minimization pro- 
cedure can stick at some of these minima before 
reaching R = 0. In that case one should think of 
the other way for solving the eq. (fl8|) . 

The ergodicity of the algorithm may also be un- 
der the question. One should check for the model 
of interest if any region of configuration space can 
be reached by the sequential updates via eq.(p"8|). 
If the algorithm is nonergodic, it can be used for 
the acceleration of other algorithms, like HMC 
and Multiboson. 

For given U, 77, to the eq.([l8|) can have many de- 
generate solutions []. This is not a problem at 

3 Such a situation takes place in SU(3) QCD, where the 



all, if one respects the symmetry relations ( |ll|jf 2j ) 
when some particular solution is being chosen. 
However, one should be careful in order not to 
violate the detailed balance and reversibility. 

In order to see, if these potential problems can 
cause any real troubles, I made some tests for the 
case of SU(3) QCD with purely fermionic term []. 
Let us however note, that these tests are prelim- 
inary, since for the reasons of simplicity and lack 
of computer time I did not include pure bosonic 
part Sb[U] into the action (this corresponds to 
/3 = 0). The more complicated simulations for 
the full QCD in physically interesting region will 
be done in the future Q . 

4. Tests for SU(3) QCD with purely 
fermionic term in the action 

The simulations were performed at 8 4 lattice 
for the partition function (Q) with 



M[U] 



(26) 



where D[U] is usual Wilson difference operator 



?(4) 



D xy =£(1- 7,0 Mx) Si + ^ y 



(27) 



and D oe means that this operator acts from even 
to odd space sites. I used even-odd precondition- 
ing to reduce the computational cost of the algo- 
rithm. By doing so one also increases the degen- 
eracy of solutions of eq. (|l8|) , which now provides 
12V* real equations for 32V^ unknown variables, 
since M[U] acts only on even space manifold 0. 

The hopping parameter k was chosen to be 
k = 0.2, which gives the plaquette value V = 
0.0089(1). The performance of algorithm was 
compared with that of usual HMC. 

The minimization of functional (19) was im- 



plemented by making the random moves in each 

expression (|is|) for unconditioned fermion matrix pro- 
vides (2 * Ndirac * N co i or * V) = 24V real equations for 
(Ng e „erator * 4V) = 32V variables. 

4 The choice of model is motivated by the popularity of 
computer simulations in QCD. 

5 One may assume, that the more degenerate the solutions 
are, the faster the configuration space is sampled, although 
this conjecture should be confirmed in the future. 
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color direction for all links lexicographically. Such 
minimization algorithm garantees reversibility of 
the procedure. 

Firstly I checked the existence of solutions of 
eq. ([l8]) at different u> for typical equilibrated con- 
figurations U. The tests were done in the inter- 
val to G [0.9, 1.1]. The good news is that for the 
considered uo one always finds some solution (i.e. 
for any 5 the minimization procedure finally gives 
R < S). No local minima of the functional ( |l9| ) 
Rmin > were observed. 

On the average solving the eq.([l8|) at e = 0.05 
with precision S = 5 * 10~ 7 required 42 min- 
imization iterations (it was checked, that im- 
proving the precision did not affect the results). 
Since the cost of one iteration is roughly equal 
to 2 * N generator = 16 multiplications by matrix 
M[{7], the total cost of finding the solution of 
eq.([l8]) was approximately 670 matrix multilica- 
tions. This has to be compared with the cost for 
generating one trajectory with HMC. At step-size 
t = 0.033 and trajectory length equal to 1 (this 
ensures 70% acceptance), one trajectory costed 
aproximately 4420 matrix multiplications i.e. 
« 6.6 times more expensive than solving eq.jl^). 

Performing the tests it was observed that our 
algorithm (let me name it Omega-algorithm in 
the following text) alone can not be used for 
sampling the partition function (|j). The prob- 
lem is that probability V[u),rj\ to accept oj < 1 
is strongly suppressed by large factor in the ex- 
ponent. One knows that the value of squared 
norm of Gaussian vector can be estimated as 
\v\ ~ Idof, where Idof is the dimensionality of 
vector space on which 77 is defined. For even-odd 
preconditioned SU(3) QCD at 8 4 lattice one has 
Idof = 24576! Therefore, using Omega-algorithm 
alone, one is restricted to choose between two un- 
favorable possibilities: 

1) Using e ~ 1/Idof- It was observed that at 
such small values of e the evolution of U fields 
becomes very slow, becouse the new configura- 
tion U always lies too close to the old one (large 
autocorrelation times); 

2) Using e 3> 1/Idof- Then u) < 1 is almost never 

6 One trajectory included 33 conjugate gradient inversions 
of matrix M[U]. Each inversion costed approximately 134 
matrix multiplications. 



accepted, and the system is gradually cooled until 
it does not move. 

In the second case the trouble appears due to 
the nonergodicity of algorithm at large e. Nev- 
ertheless, even at large e Omega-algorithm can 
be used for acceleration of other fermionic algo- 
rithms like HMC or Multiboson, since the ergod- 
icity is provided by them. 

I tested the combination of Omega-algorithm 
and HMC. After each HMC trajectory the global 
move in gauge configuration space was imple- 
mented by using Omega-algorithm at e = 0.05. 
The average plaquette value for this algorithmi- 
cal mixture was V — 0.0091(1), coinciding with 
correct one within the error bars. The efficiency 
of HMC-Omega mixture was estimated by mea- 
suring the autocorrelation times for the plaquette. 
The results were: 

Pure HMC, 1000 trajectories: r mt = 1.4(2); 
HMC+Omega, 1000 (traj + w-update): r mt = 
0.7(1). 

Note, that autocorrelation time was reduced al- 
most at no cost, becouse the global Omega- 
algorithm update is much cheaper, than HMC 
trajectory. Of course, these tests are not very 
illustrative, becouse the autocorrelation times for 
theory with pure fermionic term at this kappa 
value are rather small. The further tests in mod- 
els of practical relevance are needed. In particu- 
lar, one should find a cheap way of inclusion of 
pure bosonic action <Sf,[£/] in simulations. This 
issue is discussed in the next section. 

Let me summurize the results of the first test 
for Omega-algorithm. In its present version the 
algorithm was found to be nonergodic. Neverthe- 
less, it was demonstrated, that Omega-algorithm 
can be efficiently used for acceleration of other 
fermionic algorithms. 

Here I want to emphasize, that nonergodicity 
is not intrinsic for Omega-algorithm. It can be 
rather attributed to the unfavorable procedure 
of finding the solution of cq.(|T8|) by minimizing 
the functional (|l^), which forces us to use large 
e. One can still hope that it is possible to use 
e ~ 1/Idof, or even e = 0, and obtain nontriv- 
ial solutions U far away from U at the same time. 
This hope is based on the high degeneracy of solu- 
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tions of eq. ( fig) . The main problem is to garantee 
the reversibility of the procedure. An investiga- 
tion along these lines is in progress [H| . 

5. Inclusion of bosonic action 

Now let us consider the sampling of full par- 
tition function (||). The simplest way to include 
the contribution of purely bosonic sector is to add 
to the algorithm of the previous section the fol- 
lowing instruction: 

• accept the new configuration U with the 
probability 

P*° s c on = mm (1, e Sb[u] ~ Sb[G] ) (28) 

The reader can easily check that one gets a valid 
algorithm for sampling the partition function (|J) ■ 
However, the new configuration U, obtained by 
solving the eq. (filf) , differs globally from the old 
one, so the acceptance ( p8| ) can be very small. Of 
course, one can try to make the new configuration 
lying close to the old one to ensure the reasonable 
acceptance, but by doing so one lose one of the 
main advantages of the new algorithm - the finite 
step-size of update. 

A more radical way is to rewrite the partition 
function (||) in the form: 

^./ W (#RW|)i7 (29) 

where operator B satisfies the following identity: 

Det B[U] = e- Sb[u]/2 Det M[U] (30) 

After that one can use the algorithm of the pre- 
vious section with the equation 

B[U]r) = uj B[U]r] (31) 

instead of eq. ( fl8| ) . There is a large degree of am- 
biguity in defining the operator B. For example, 
one can choose B[U] = e~ Sb ^/ 2N M\U], where N 
is the size of matrix M. A more rational way is, 
probably, to insert the bosonic contribution into 
the operator M locally. One can represent the 
bosonic action as a sum over local contributions: 

1 N 

-S b [[/]=5>«[[/] (32) 
i=l 



(some of S¥' can be equal to zero), and then de- 
fine the matrix B as follows: 

B i:j = e- s " Mij l<i,j<N (33) 

(no summation over i is assumed in ([53])). In any 
case, one should try to choose the most conve- 
nient version of operator B in order to simplify 
the procedure of solving the eq.(|3l|). 

6. Discussion. 

The new method, proposed in this paper, can 
provide a cheap simulation algorithm for theories 
with dynamical fermions. The algorithm is exact 
and has a finite step-size of update. 

Although in its present version the algorithm 
was found to be nonergodic, it can be used in 
combination with other ergodic algorithms like 
HMC and Multiboson. One may expect a good 
performance for such algorithmical mixture when 
the lattice volume increases, since computational 
cost for our new algorithm grows only lineary 
with the volume of system. 

The efficiency of algorithm can be strongly af- 
fected by the clever choice of procedure of solv- 
ing the eq.(p^|). Due to the high degeneracy of 
solutions one may hope to find the procedure, 
which makes the algorithm ergodic, and samples 
the configuration space fast enough. Also one can 
speculate on possibility of finding some analytic 
solutions for the eq.JTsI). If this is feasible at all, 
one can obtain a very cheap fermionic algorithm, 
comparable in cost with the algorithms for purely 
bosonic theories. 
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